Josephson vortices and the Meissner effect in stacked junctions and layered 
superconductors: Exact analytical results 
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We present an exact mathematical description of Josephson vortices and of the Meissner effect in 
^ ' periodic thin-layer superconductor/insulator structures with an arbitrary number of identical junc- 

tions TV — 1 (2 < TV < oo, where TV is the number of superconducting layers) in terms of localized 
■ solutions to a system of differential equations for phase differences. We establish a general criterion 

of the existence of localized solutions. We show that Meissner solutions are characterized by several 
Josephson lengths Xji lengths for even TV, and ^fci lengths for odd TV). We derive an exact 
expression for the superheating field of the Meissner state, H s , as an explicit function of TV. For 
Josephson vortices, we find two basically different types of topological solutions: "vortex-plane" 
solutions and incoherent vortex solutions. Thermodynamically stable "vortex-plane" solutions rep- 
resent a chain of TV — 1 vortices (one vortex per each insulating layer). They are characterized 
^ , by the same set of Xji as the Meissner solutions. We obtain exact analytical expressions for their 

S_h ■ self-energy and for the lower critical field H c \. Incoherent vortex solutions comprise solutions with 

Oh' k < TV — 1 vortices and different vortex-antivortex configurations. In contrast to the "vortex-plane" 

solutions, they prove to be thermodynamically unstable, and their spatial dependence is character- 
ized, in general, by TV — 1 length scales. As an illustration, we analyze 1-4-Josephson-junction stacks 
and investigate a transition to the layered superconductor limit (TV — » oo). 
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O ; I. INTRODUCTION 

We present a rigorous mathematical examination of the problem of Josephson vortices and of the Meissner effect 
in thin-layer Josephson-junction stacks and layered superconductors, with a static external magnetic field H applied 
parallel to the layers (along the z axis, see Fig. 1.) We consider periodic systems composed of an arbitrary number 
j^f~\ ' TV — 1 of identical superconductor/insulator (S/I) junctions (2 < TV < oo, where TV is the number of S-layers, with 
x being the layering axis). Our starting point is the microscopic Gibbs free-energy functional derived in Ref. [El. 
<r— ( ■ Mathematical structure of this functional is analogous to that of the phenomenological Lawrence-Doniach model. B] 
Thus, the treatment of our paper fully applies to the latter model as well. 

Mathematically, both Josephson vortices in moderate fields |H| and the inhomogeneous Meissner state are de- 
scribed by solutions of a system of nonlinear second-order differential equations for phase differences, (f> n , with 
square-integrable first-order derivatives. (For brevity, we call here such solutions "localized"). Our approach 
is substantially based on the observation of a nontrivial property of the differential equations for <j) n : We show that 
the problem of finding localized solutions for <j) n can be reduced to solving a standard initial value problem. Using this 
key mathematical result, we establish an exact criterion of the existence of localized solutions. The existence criterion, 
O ■ in turn, allows us to obtain a complete classification of physical localized solutions. We have found three types of such 
solutions: Meissner solutions, topological "vortex-plane" solutions, and topological incoherent vortex solutions. 

Meissner solutions are localized near the side boundaries y = —L and y = L. In contrast to the well-known single- 
junction case, HH the Meissner solutions in stacks with TV > 4 turn out to be characterized by several Josephson 
lengths Xji lengths for even TV, and lengths for odd TV). The Meissner solutions persist up to a certain 

superheating field of the Meissner phase, H s . We derive an exact expression for H s as an explicit function of TV. We 
show that the field H s simultaneously determines the penetration field for "vortex planes". (See below.) 

Thermodynamically stable " vortex-plane" solutions represent a chain of TV — 1 Josephson vortices (one vortex per 
each I-layer), positioned in the symmetry plane y — 0. These solutions are uniquely determined by the vortex- 
penetration conditions at |H| = H s . They are characterized by the same set of Aji as the Meissner solutions. Such 
solutions were previously obtained for infinite (TV = oo) layered superconductors . [p],p| Under the name of the "coherent 
mode" or the "in-phase mode" they are well-known in double-junction stacks. (4|-|7|] For 4 < TV < oo, the existence of 
the " coherent mode" was predicted in Ref. j| . (The authors of Ref. Q specially emphasized the importance of this 
mode for practical applications.) Besides giving a proof of the existence and stability of the "vortex planes" in the 
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general case 2 < A < oo, we derive exact analytical expressions for their self-energy E v and for the lower critical field 

Incoherent vortex solutions comprise single- vortex solutions, vortex solutions with 2 < k < N — 1 vortices in the 
plane y — 0, as well as different vortex-antivortex configurations. All such solutions satisfy the existence criterion. 
However, in contrast to the "vortex-plane" solutions, they prove to be thermodynamically unstable and do not meet 
the vortex-penetration conditions at any |H| ^ 0. It should be noted that the single- vortex solutions obtained 
in this paper have no resemblance to hypothetical Abrikosov-type vortices, introduced without proper mathematical 
justification in some previous publications. |lO| , |ll"| ] Besides being thermodynamically unstable, the actual single- vortex 
solutions are not uniquely determined by asymptotic boundary conditions. They are accompanied by singular phase- 
difference distribution in all A — I junctions, and their spatial dependence is characterized, in general, by A — 1 length 
scales. 

Section II of the paper is devoted to exact mathematical formulation of the problem. In section III, we derive all 
major physical and mathematical results sketched above. The general consideration of this section is illustrated by 
several concrete examples in section IV. In particular, we analyze 1-4-junction stacks and investigate a transition to 
the layered-superconductor limit (A — > oo). The obtained results are discussed in section V. Appendices A-C contain 
some additional mathematics, relevant to the subject of our study. 



II. FORMULATION OF THE PROBLEM 



We begin by writing down the microscopic Gibbs free-energy functional |ffl] of a periodic structure consisting of 
alternating A superconducting (S) and A — 1 insulating (I) layers (2 < A < oo): 
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<Pn(y) = <Pn(v) ~ <fn-l{y)- 

Here ft = c = 1; a is the S-layer thickness; p is the period, and W z is the length of the structure in the z direction 
(W z — > oo); the length of the structure in the y direction is W y — 2L; f n (y) [0 < f n (y) < 1] and <p n {y) are ; respectively, 
the reduced modulus and the phase of the pair potential A n (y) in the nth superconducting layer: 
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A n (y) = A(T)/„(y)exp<p„(y), 

with A(T) being the microscopic gap at temperature T; £o is the BCS coherence length; £(T) and X(T) are, re- 
spectively, the Ginzburg-Landau (GL) coherence length and the penetration depth; Z^cos^) is the incidence-angle- 
dependcnt tunneling probability of the I-layer between two successive S-layers; H C (T) is the thermodynamic critical 
field; A = (A x ,A y ,0) is the vector potential. The local magnetic field H(x,y) = [0, 0, H(x, y)] obeys the Maxwell 
equation 

dA y (x,y) dA x (x,y) 



dx dy 



H{x,y) 
with boundary conditions 

H(Q,y)=H((N-l)p,y)=H, y&[-L,L], 

H(x,±L)=H, xe[0,(N-l)p], 

where H is a static external magnetic field applied along the z axis. (See Fig. 1.) The sum of the first three phase- 
and field-independent terms on the right-hand side of ([!]) represents the condensation energy. The fourth term is 
the kinetic energy of the intralayer currents. The last two terms are the Josephson energy and the field energy, 
respectively. 



Expression (Q) is valid under the conditions 



T -^r L « i, (2) 



£o « a, (3) 
a«min{C(T),A(T),a- 1 e }, (4) 

a<Cp. (5) 

Conditions (|^) (T c q is the critical temperature of an isolated S-layer) and (|3|) ensure the applicability of the GL-type 
expansion within each S-layer. Condition (^) corresponds to the thin S-layer limit, whereas condition ([5]) is employed 
here for the sake of mathematical simplicity only. Being a first-order expansion in a/p, equation ([!]) applies in fields 
\H\ <S H C 2, where H C 2 is the upper critical field. 

Mathematical treatment of the functionals of the type (0) is described in full detail in Ref. ||), section III: One 
minimizes (|l|) with respect to /„ and A x , A y , imposes the gauge A x — and eliminates A y by integration. The result 
is a closed, complete set of coupled nonlinear mean-field equations for the reduced modulus of the pair potential /„ 
and the phase differences </>„, together with relations for all physical quantities of interest. To simplify mathematical 
analysis of the mean-field equations, we introduce dimensionless units by 

x 

> x, 

P 



Voo 



H 



where the quantities on the left-hand side are dimensional, with Xjoo = (8nejop) being the Josephson penetration 
depth (jo is the density of the Josephson current in a single junction with thick electrodes) and H soo — (epXjoc) 
being the superheating (penetration) field of the infinite layered superconductor. 0,^] In our dimensionless units, 
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for example, the flux quantum is $o = and the lower critical field of the infinite layered superconductor JllEI is 
u — 1 

In the dimensionless form, the mean-field equations for /„ and 4> n read 



foiv) ~ foiv) = r(T) 



e 2 d 2 f (y) 2[H-H 1 {y)f 1 



fn(y) - f n {y) = r(T) 



e 2 rf 2 / n (y) 2[g n (y)-g n+1 (y)r 
2 dy 2 + e2/ 3( y) 
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1 < 71 < AT - 2, 
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' [H n+1 (y) - H n (y)} - [H n (y) - ff„-i(y)] - e 2 H n (y) 



my) 



fLM 



where 



d(j) n 
dy 



-t^M, Kn<N-l, 
2 dy - - 

H (y) = H N (y) = H, 
(±L) = 2H, l<n<N-l, 



< 1, 



A 



and the local magnetic field in the nth insulating layer (n — 1 < x < n) is given by 

v 

H n {y) = - / duf n (u)f n -i(u) sin (f> n (u) +H 



(8) 
(9) 
(10) 

(11) 



1 

- J duf n {u)fn-x{u) sin </>„( 



(12) 



Note that although in most physical situations e <C 1, for the sake of mathematical generality, in this paper, e is 
supposed to satisfy only the weak inequality (|Tl|). Because of the obvious property f n (y) — fn{~y)j equation ( O ) 
implies that the phase differences 4> n meet the condition 
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<j> n (y) = -M-V) + 0mod2^. (13) 
The dimensionless Gibbs free energy fl(H), normalized via the relation 

H*(T)a\ Joo W, n[Hh 
in terms of the mean-field quantities f n , <j) n and H n (y) has the form 

L 



N-l J 

n w = E / d y 



-/^) + k(y) + r(T)e 



2 J ""' y 2 V rfy 



2ri 1 ' [fln+i(y) - ^(y)] 2 + 2r(T) - tf] 2 



^ 2 /„ 2 (y) 



+ "V E 7 d » + fniv) - 2fn(y) f n-l(y) COS <t> n (y)] 



(14) 



Here, the two terms in the second line on the right-hand side are the kinetic energy of the intralayer currents and the 
field energy, respectively. The intralayer current in the nth S-layer J n (y) (normalized to H soo ) and the density of the 
Josephson current between the nth and the (n — l)th S-layers jn,n-i(y) (normalized to jo) are given by 

J n (y) = ^-[H n (y)-H n+1 (y)}, < n < N - 1, (15) 

and 

in,n-i(y)=2^M = / n (y)/„_!(y) sin <t> n (y), 1 < n < N - 1, (16) 

respectively. Relations (p|)- (|l"6| ) provide a complete, self-consistent description of the thin-layer periodic S/I structure 
in the temperature range (0) and in fields |if| <C H C 2. 

In this paper, we will be interested in physical solutions for <f> n with a square-integrable first-order derivative, 
localized within a spatial range of order unity. (For brevity, we call such solutions "localized".) Therefore we assume 
the condition 

£>1. (17) 

Moreover, we assume that the temperature range satisfies the condition of the weak-coupling limit 

r(T) < 1. (18) 

One can obtain a perturbative solution for /„ and 4> n up to any desired order in r(T), starting from the zero-order 
solution to (||), (|7|), 

fn = 1, (19) 

and the zero-order equations for <f) n , 

H n+1 (y)-(2 + e 2 )H n (y)+H n _ 1 (y) = ~^l, 1 < n < N - 1, (20) 

where H n {y) are given by (^) with f n = 1 and satisfy the boundary conditions (^J). For most applications, it is 
sufficient to consider expressions for physical quantities only in leading order in r(T). Thus, for example, substituting 
( |l9| ) and the solution of ( pfl| ) into (|l4|) immediately yields a first-order expansion for the Gibbs free energy, because 
first-order corrections to the condensation-energy term cancel out. 
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A detailed mathematical analysis of Eqs. (^) is the subject of section III. Here we point out that these equations 

can be transformed into a very useful for application form by solving for H n (y) (see Appendix A for mathematical 
details): 

H n (y) = h n (y) + H n , (21) 

M») = y£G(»,mA^, (22) 

m=l " 

^" - ,,-N _ n N ' W 



where G(n,m) are given by (A9), and /x is given by ( A5). By (AS), and (JlCj), ( A15 ), expression (|2l| ) explicitly 



satisfies boundary conditions (|9|) and H n (±L) = H. Moreover, the ^-independent quantities H n in (|21j) have clear 
physical meaning: Being solutions of ( |20| ) with d ^2y V ^ = 0, they describe distribution of the local magnetic field 
within I-layers in the homogeneous Meissner state (see section III of Ref. ||). Also note that H n = i?jv-n> which 
is a reflection of the symmetry of the problem. [By comparison, in an infinite layered superconductor H n = 0, and 

N—l +oo 



y~] G(n, m) . . . — > ^2 G co (n, m) where Goo(n, m) are defined by (A15).] 

m—l m— — oo 

In addition, we point out that equations of the popular phenomenological Lawrence-Doniach model [Q also can be 
reduced to the dimensionless form (|6|)-(|l0|), jT^)-([l^), with r(T) being a phenomenological parameter, e = and 

Q(H) normalized to gc(T) ^ Jcoff ' . (See Ref. || for more details.) Thus, all the consideration of this paper fully 
applies to the Lawrence-Doniach model as well. 

III. MAJOR RESULTS 
A. The criterion of the existence of localized solutions 



By differentiation with respect to y, integrodifferential equations ( |20| ) reduce to a system of A — 1 ordinary nonlinear 
second-order differential equations 

= 1 [(2 + e 2 ) sin MV) ~ MV)} , 
d ^ = ^ [(2 + e 2 ) sin (j) n {y) - sin0 n+1 (t/) - sin^_i(y)] , 2 < n < N - 2, 

d2<t)N d ^ iv) = ^ [(2 + ^ ^MMv) - smMMv)} (24) 

with boundary conditions (|l0|). 

Consider Eqs. ( |24| ) on the whole axis — oo < y < +00. Two simple properties of ( g4j ) are quite obvious: If <p n (y) 
(1 < n < N — 1) is a solution, the functions 0„(y) given by 

0n(y) = My) + 2ttA: (fc is an integer), (25) 

and 

4>n{,y) = 4>n{y + c) (c is an arbitrary constant) (26) 

are also solutions. [The latter is a result of the fact that y does not enter explicitly the right-hand side of (pjj).] Our 
conclusions about the existence of localized solutions to ( |24| ) will be substantially based on another key property, 
which we formulate as a lemma: 
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Lemma. Consider an arbitrary interval I = \L\,L^ and yo € /. The initial value problem for Eqs. ( pi] ) with 
arbitrary initial conditions <p n {yo) = -^(yo) = Pn has a unique solution in the whole interval /. This solution 
has continuous derivatives with respect to y of arbitrary order and continuously depends on the initial data. (For the 
proof of the Lemma, see Appendix B.) 

It is worth noting that the existence and uniqueness of a smooth solution to the initial value problem in the whole 
interval / is rather nontrivial for nonlinear differential equations: For such equations, theorems of existence and 
uniqueness are usually valid only locally, in the neighborhood of initial data. Jl^] In our case, global character of the 
solution and its infinite differentiability are ensured by the fact that <fi n enter the right-hand side of Eqs. (|I|) only as 
arguments of the sine. Note that because of the arbitrariness of the interval /, the solution can be uniquely continued 
onto the whole axis —oo<y< +00. |Q Now we will show that the problem of finding localized solutions to (|24|) can 
be reduced to the standard initial value problem. 

Differentiating (|2l]) with respect to y yields 

sin 0„ (y) = e 2 £ G(n, m) . (27) 

m=l " 



Multiplying ( p7| ) by summing over the layer index n with the use of ( A10| ) and performing integration, we 

arrive at the first integral of Eqs. (p3): 



N-l 2 w-i N-l 



n ST x t \ e r>( ^d(j) n {y) d(j) m {y) 
C-^ C0S( l> n { y ) = -^^G(n,m)— — , (28) 

n=l n=l m=l y y 

where C is the constant of integration. Now let us choose an arbitrary point yo € [-L, L], where L is sufficiently large 
[see the condition (fL7|)]. We are looking for localized solutions of Eqs. ( pi| ) that in the region 

A roax < \y - y \ , (29) 

where A max is determined by the maximum positive eigenvalue of the symmetric matrix G(n, m) (see Appendix A), 
satisfy the asymptotic conditions 

<f>n(y) =0mod27r + o(l), (30) 

— r =o{l) (31) 

for any 1 < n < N — 1. By inserting ( jfidj ) and (^) into d2S|), we establish the value of the constant of integration: 
C = N — 1. Thus Eq. (H) becomes 

n=l n=l m=l a a 



Substituting initial values 4> n {yo) — &m "r^iVo) — Pn mto (p2|), we obtain the general criterion of the existence of 



dy 

localized solutions to (]24|): 



N-l 2 N-l N-l 



^ Sin2 T = T ^ ^ G(n,m)0 n p m . (33) 

n— 1 n—1 m— 1 

Indeed, the Lemma guarantees the existence, uniqueness and differentiability of a solution for arbitrary a n and /? n in 
the whole interval [— L, L\. Owing to our special choice of the constant of integration in (|3^), the solution determined 
by a n and j3 n obeying ( |33| ) will necessarily satisfy asymptotic conditions (|o|), ( |3l| ) in the region (^9|). Moreover, this 
solution will automatically satisfy the conditions 

(34) 

in the region (E9) for any I < n < N — 1 and 2 < k, which can be verified by repeated differentiation of d2^) and 



application of (p0[), (31). Below, we apply the criterion (p3[) to concrete physical situations. 
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B. Meissner solutions. The superheating (penetration) field H s 

Let H > for definiteness. Consider an intermediate length scale R such that 

!<*<§. 

(See Fig. 1.) The Meissner boundary value problem is specified by the boundary conditions ( |l0| ) and 

<P n (±(L-R))=o(l), ^(±(L-R))=o(l) 

dy 

for any 1 < n < N — 1. For — L < y < —L + R, we write using ([l2|): 

-ffn(y) = ^ y cftt sin <^„(u) + - J du sin <j) n (u) + H. 

-L+R -L 

According to (|2|), we make the identification 

y 

My) = 2 



(35) 



(36) 



(37) 



eh/ sin (f> n (u), 



(38) 



-L+i?. 



Hn — — 

2 



J dy sin 0„ (y) + H, 



(39) 



where /i„(?/) stands for the field penetrating through the y = —L interface. Using the fact that H„ < H, and 
H n = Hn-ui we establish the following important properties: 



7T < My) < o, 



(40) 



My) = <pN-n(y), K{y) = h N My), H n {y) = h n Mv)- ( 41 ) 

Relations ([ill ) are a result of the symmetry of the problem. They imply that the number of independent equations 
describing the Meissner solution is ^ f° r even N and for odd N. The solution for <p n (y) in the region L—R < y < L 
can be obtained from the solution in the region — L < y < —L + R using the property 

My) = -M-y), (42) 



resulting from the general relation (13). In the region —R < y < R, <j> n = 0, and we have 

H n (y)=H n . 



(43) 

Now we apply relation (|3^) for yo = —L. In view of the boundary conditions ^-(—L) = (3 n = 2H [sec (ffoh], we 



get 



N-l 



N 



1 \ . 9 a 

> sin — 

- 1 ^ 5 



IP 

m ' 



(44) 



where 



H s = 



1 - 



V 1 + T-^J (l-M"- 1 ) 
e(iV- 1) (l + ^-i) 



(45) 
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Physical interpretation of the quantity H s is straightforward. The maximum value of the left-hand side of ( f44| ) is 
unity, which corresponds to H = H s on the right-hand side. Hence H s is the maximum external field for which a 
Meissner solution is still possible, i.e. the superheating field of the Meissner state. ^,|| Moreover, the maximum of 



the left-hand side is achieved when all a n = cj> n (— L) = — it, which, by (42), gives 4> n {L) — (f> n (—L) — 2-k for all n. The 
total flux of the penetrating field at H = H s 



N—l ~ L + R N-l u 

$ = E J d v h n(y) + E J d y h n(v) 



= ir(N-l) 



2J1 



/' 



JV-l 



e(N-l) l + 



n 



is exactly equal to the total flux carried by a "vortex plane", i.e. a chain of N — 1 vortices positioned in the plane 
y = 0. [See Eq. ( |52|) below.] These conditions correspond [^]|l] to simultaneous and coherent penetration of Josephson 
vortices into all the junctions. Therefore H s can be also regarded as the penetration field for a "vortex plane". Note 
that H s for N < oo is always higher than the penetration field of an infinite layered superconductor H soo — 1 . Jl],|| 

Thus, the Meissner boundary value problem, Eqs. (|l(]) and (|36|), for H = H s reduces to the initial value problem 
with a„ = cj) n {—L) = — 7r and f3 n = d< ^ n jy L ^ = 2H S (1 < n < N — 1). Let ip^'iy) be the corresponding unique 
solution. By continuous dependence of a solution of the initial value problem on initial data (see the Lemma), we 
can argue that the existence of the unique Meissner solution <j>^ a (y) at H = H s guarantees the existence of a unique 
Meissner solution at any H < H s with a certain set of a n (— n < a n < 0) satisfying (f44|). In other words, we have 
proved that the Meissner boundary value problem, Eqs. ([l0|) and (^), has a unique solution for any H < H s . 

Finally, we want to emphasize one more important property of Meissner solutions in (N — Injunction periodic 
S/I structures that was not realized in previous publications. Owing to the fact that these solutions satisfy a 
system of coupled second-order differential equations (-j equations for even N, and equations for odd iV), they 
are necessarily characterized by several Josephson lengths Xji (^ lengths for even N, and ^ lengths for odd N), 
in contrast to a single junction with only one Aj. All Aj,; are determined by positive eigenvalues of the symmetric 
matrix G(n, to). (See Appendix A.) 



C. Vortex-plane solutions. The lower critical field H c i 

By a vortex-plane solution we understand a chain of N— 1 Josephson vortices (one vortex per each I-layer) positioned 
at y = 0. Analogously, an antivortex-plane solution is a chain of N ~ 1 Josephson antivortices (one antivortex per 
each I-layer). Such solutions are characterized by the symmetry 

Mv) = ± 2 *-M-v) ( 46 ) 

[see (fl3|)] and asymptotic boundary conditions [|l3| 

(j> n {-R) = o(l), <j) n {R) = ±2ir + o(l), (47) 



for all 1 < n < N — 1 and any k > 1. [The "plus" sign in ( |4q ) and (47) corresponds to a vortex plane in fields H > 0, 
whereas the "minus" sign corresponds to an antivortex plane in fields H < 0.] The existence and uniqueness of these 
solutions follows immediately from the Lemma, the criterion ( |33| ) and the results of the previous subsection. Indeed, 
by (f4(|), a vortex (antivortex) plane satisfies the conditions a„ =</>„(0) = ±7r for all n. Under these conditions, the left- 
hand side of ( |33| ) reaches its maximum. As shown in the previous subsection, the maximum condition corresponds 
to the unique choice j3 n = ^p-(Q) = ±2iJ s on the right-hand side of (|33|), consistent with the vortex-penetration 
conditions (see the next paragraph). Hence the initial values 
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,(0) = ±7T, 



dy 



±2H„ 



1 < n < N - 1. 



(49) 



determine a unique localized solution in the region —R<y<R that automatically meets the asymptotic boundary 
conditions @, ©. 

Thermodynamic stability of vortex-plane (antivortex-plane) solutions is ensured by the fact that they satisfy the 
vortex-penetration conditions for H = ±H S . Let H > for dcfiniteness. (In what follows, we consider only vortex 
planes. The discussion of antivortex planes is quite analogous.) A vortex-plane solution can be constructed from 
the Meissner solution cf>^ s in the region — L < y < —L + R, discussed in the previous subsection. Indeed, using the 
properties (p5~|), (p6|), we obtain a solution 

&(y) ^ <P" s (y- l) + 2tt 

in the region < y < R that satisfies the initial conditions ^„(0) = tt, ^ L (0) = 2H S . This solution can be continued 
p2[ into the region — R < y < 0. By the uniqueness of a solution to the initial value problem, the obtained solution 
coincides with the vortex-plane solution in the interval —R < y < R. 

The local magnetic field in the presence of a vortex plane is given by general relations (^l|) and (|37|)-(p9|), where 
h n (y) [with y £ (~R, R)] should be interpreted as the magnetic field induced by the vortex plane itself. Owing to the 
vortex-plane initial conditions (3 n =-^-(0) = 2H S , 



MO) = H s 



-N+n 



N-nl 



-N 



,N 



(50) 



A 4 - V> 

(See Fig. 2.) The symmetry relations (Ml|) apply to vortex-plane solutions too. Thus, the number of independent 
equations describing a vortex plane is for even N, and for odd N, as in the case of Meissner solutions. 

Consequently, spatial dependence of vortex-plane solutions is characterized by the same set of Xji as that of the 
Meissner solutions. 

The flux <E> n through the nth I-layer can be found using (p|), (E^) and (|A13|): 



R 



.-N+n 



®n = J dyh n (y) = tt 

-R 

The total flux carried by a vortex plane is 

N-l 
n=l 



1 - 



M A 4 



N—n "I 



-N 



(51) 



iV 



e(iV-l) 



1 



(52) 



Note that in contrast to Josephson junctions with thick electrodes [g[ and infinite layered superconductors, |1|J^] the 
flux carried by a Josephson vortex in a finite thin-layer S /I structure is not quantized and is always smaller than the 
flux quantum $o = tt. (This fact has been already pointed out in Ref. 0.) 

To determine the thermodynamic lower critical field H c \ at which the vortex-plane solutions become energetically 
favorable, we must calculate the difference between the Gibbs free energy in the presence of a single vortex plane, 
ily(H), and the Gibbs free energy of the homogeneous Meissner state, £Im(H) [the sum of phase-independent terms 
in Q]. Substituting @-@ into Q and using (|A7|), @, in first order in r(T) < 1, we obtain: 



n v (H) - n M (H) 



r(T) 



N-l N-l 



^2 X! G ( n i m ) / d y 

— i ™ i J 



n—l m— 1 



d<j) n (y) d<fi m (y) 
dy dy 



A<S>H 



(53) 



where the total flux <!> is given by (|52j). The first term on the right-hand side of ([53|) should be interpreted as the 
self-energy of the vortex plane: 



N-l N-l 



1 ^ i J 



dy 



n—l m— 1 



d<t> n (y) d<j) m (y) 
dy dy 
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J N-l 

4r(T) dy^2 sin 



><Pn(y) 



r(T) 



N-l 

E 

n=l 



dy 



N-l 



— G(n,m) 



d(j>n{y) d(/> m {y) 
dy dy 



+ 1 - cos(fi n (y) 



(54) 



Note that E v is exactly twice the energy of the Josephson currents in (|14j). Besides, formulas (|53|) and (]54|), 
with corresponding reinterpretation of fl v (H), $ and £J„, also hold for incoherent vortex solutions considered in 
the next subsection. They also apply to an infinite layered superconductor, taking account of the substitution 

N-l +00 

G(n,m) . . . — ► Goo(n,m) . . ., where Goo(n,m) are defined by (A15). In this latter case, the self-energy 

m,n— 1 m,n— — 00 

should be additionally minimized with respect to the phases ip n , which immediately yields the exact solution [QJ|] 
with 4> n (y) = (f>n+i(y) = 4>{y) an d A joo = 1. [In the case N < oo, the minimization with respect to ip n is not allowed 
by the boundary conditions (||). 
From (p3|), we get: 



J5« 



4r(T)$' 



(55) 



For thermodynamically stable solutions, we must necessarily have H c \ < H s . It is straightforward to verify that this 
condition is met by the vortex-plane solutions. Using (p7|), (Efj) (with the "plus" sign), the initial values f3 n = 2H S 
and integrating by parts, we convert E v into the form 



E v = 2r(T) 



jv-i „ 

Y] / dy(j) n (y) sin cj) n (y) - 2H S § 



The first term on the right-hand side of (pfl ) is positive, because in the region < y < R all 
7T < <f) n < 27T- By the use of (|2^ ) and ( |A14| ), we obtain the following strict inequalities: 



(56) 

satisfy the relation 



JV-l 



n=l 



dy0„(y)sin^ n (y) < 4i7 s $, 



< E v < Ar{T)H s §. 

Hence, 

< Ha. < H s , 

as anticipated. Note that in all special cases admitting exact analytical solutions (N — oo, jl],[5| and N = 2,3, see 
section IV), H c \ = —H s . Finally, we want to point out that in contrast to vortex-plane solutions in infinite layered 
superconductors, where H n (y) — H n+ i(y) — H{y) and J n = for all n, in finite structures the intralayer currents 
J n , in general, are not equal to zero, as can be easily seen from (|l5|). Only for even number of junctions (N is odd), 
in the central S-layer J «-i = 0, by the symmetry (plj) . 

D. Single- vortex solutions and other localized incoherent vortex solutions 

A single Josephson vortex positioned in the Zth I-layer at y = obeys symmetry relations 

(fii(y) = 2tt - <fii{-y); <fi n {y) = -(fin(-y), n^l, (57) 
[see (|||)] and asymptotic boundary conditions [|l3| 

fai-R) = o(l), MR) = 2tt + (1), (58) 
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<f> n (±R) =0(1), n^l, 



(59) 



rf fc <M±i?) 

dy k 



o(l), for all 1 < n < N — 1, and any k > 1. 



(60) 



Moreover. 



dh n ( y ) 

dy 



> in the region — i? < y < 0, and dh " l - y ' < in the region < y < R. Hence, <j) n must satisfy the 



dy 



relations 

< <p n (y) < ir, y S (-i?,0); -tt < 0„(y) < 0, y S (0, R) , for n ^ I, 
[see ( p8| ) for y £ {—R, R)] and the initial conditions 

ai = 4>i(Q) = n; a n = 0„(O) = 0, n^l, 



(61) 



(62) 



_ d0i(O) #»(0) ^ n , , 

Pi = — j > 0; p n = — < 0, fif i. 

dy dy 

A necessary condition of the existence of such solutions is provided by the general criterion 

9 N-l N-l 



J2 G{n,m)l3 n l3 m = 1. 



(63) 

|) and has the form 

(64) 



n—1 m— 1 



Relation ( |64| ) imposes only one constraint on N— 1 quantities (3 n and can be satisfied by different sets of /?„. Therefore, 
in contrast to the vortex-plane problem, a solution to the single-vortex problem for any given vortex position I is not 
unique. The determination of the optimum set of (3 n may require an additional examination of the vortex self-energy 
(see section III.B). 

Besides failing to meet the requirement of uniqueness, single- vortex solutions, in general, break the symmetry (|4l|), 
inherent to the original integrodifferential equations (^), (|2^) minimizing the Gibbs free-energy functional ([!]). (The 
only exclusion is a stack with an odd number of junctions N — 1 and I = Even more important is the fact that 
single-vortex configurations do not satisfy the vortex-penetration conditions for any H > 0. Indeed, starting from a 
single- vortex solution in the region —R<y<R,we can construct a solution in the region — L < y < —L + R satisfying 
the Meissner boundary conditions (36). (Compare with the reverse procedure of the construction of a vortex-plane 
solution from the Meissner solution (/>n s > described in the previous subsection.) However, the solution thus obtained 
will not represent any physics, because it cannot meet the physical boundary conditions d( ^" c |~ L - ) = j3 n — 2H for any 
H > 0: According to (|63|), fli and [3 n with n 7^ / have different signs. Physically, this means that isolated vortices 
cannot penetrate the periodic S/I structure at any static H > 0, and the penetration field for isolated vortices cannot 
be defined. On the basis of these observations, we conclude that, in contrast to vortex-plane solutions, single-vortex 
solutions are thermodynamically unstable and do not represent any solutions to (|^), (^0|) for H > 0. This situation has 
a simple mathematical explanation. In contrast to differential equations (|24|), integrodifferential equations (||), j20| ) 
explicitly contain the external magnetic field H . [See the explicit expressions for H n (y), Eqs. (jl^).] In the case of Eqs. 
(p4|), the external field H enters only via the boundary conditions (|l^) at y = ±L, whereas single- vortex configurations 
are required to satisfy asymptotic boundary conditions (p8|)-(|60|) at y = ±R. It is therefore not surprising that Eqs. 
(p4|), restricted to the interval (— R, R) C [— L, L], may possess redundant solutions that do not satisfy Eqs. (||), (p?]). 
[Note that the exact equations (||), valid for arbitrary r(T), do not, in general, reduce to any differential equations, if 
fn{y) ^const.] 

The flux through the nth I-layer due to the vortex in the Ith I-layer can be found using (E2n and (pq), (J59|): 



4>, 



dyh n {y) 



-R 



2^ I 



\ n -i\ 



The total flux carried by the vortex in the Zth layer is 



l-N 



,N-l\ 



1 - 



N-r 



li- 1 + ii- N+l n N ~ l 



(65) 



(66) 
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For N < oo, the total flux <I> is not quantized and is less than the flux quantum <f>o = tt- (See the previous subsection.) 
Note that the total flux carried by a single vortex in the Zth I-layer is exactly equal to the flux through the Ith I-layer 
in the case of a vortex plane, given by Eq. (51) with n = I. 

The consideration of other localized incoherent vortex solutions to ( p4| ) (i.e., solutions with 2 < k < N — 1 vortices 
in the plane y = 0, and vortex-antivortex configurations) can be done along the same lines. All these solutions are 
thcrmodynamically unstable too. We want to underline that our general conclusions about the form of static single- 
vortex configurations completely agree with the results of numerical calculations of Ref. |Q. (See, in particular, Fig. 
5 therein.) On the other hand, these single-vortex configurations have no resemblance to hypothetical Abrikosov- 
type vortices introduced without appropriate mathematical justification in Refs. jn^fOJ . Let alone thermodynamic 
instability, the actual single-vortex solutions to (|24|) are accompanied by singular phase-difference distribution in all 
N — 1 junctions, satisfying (|6l[)-(p3|) and the existence condition (p4[). Their spatial dependence is characterized, in 
general, by N — 1 length scales, which is inherent to localized solutions of a system of coupled second-order differential 
equations. These intrinsic features of topological solutions in discrete periodic S/I structures cannot be reproduced 
by any imitation of Abrikosov vortices, typical of continuum type-II superconductors. It is also worth reminding that 
Abrikosov vortices in the London approximation are described by linear partial differential equations, |H| whereas 
the ordinary differential equations (24|) are essentially nonlinear. The principle of superposition of solutions is not 
valid for nonlinear equations. Unfortunately, this basic point is sometimes disregarded in literature. |l4| 



IV. PARTICULAR EXAMPLES 
A. A single thin-layer junction (N = 2) 



In this simplest case, the only nonzero element of the matrix G(n,m), given by (A9), is 



G(M) = 2T^- (67) 



By (27), a single phase difference 4>i(y) = 4>(y) satisfies the usual static sine-Gordon equation 

d 2 Hy) 1 



= —75- sin 



dy 2 A 2 
with the Josephson length @] 



V2 + e 



<Kv), (68) 



(69) 



2 



Note that Aj, given by (j69|), for e < 1 is much smaller than the Josephson length of a single junction with thick 
electrodes, which in our dimensionless units is A,/o = \/3x" @ From (^3|) and (|45|), we get the local field in the 
homogeneous Meissner state 



V2 



e 



2 



and the superheating (penetration) field 



respectively. For e -C 1, the superheating (penetration) field H s , given by (|71[), is much higher than the corresponding 
field j^J^] of a single junction with thick electrodes H s q — Xjq. 



Hs = Xj L = (71) 



1. The Meissner solution 



For the fields < H < H s = y ' 2 ^ e2 , the Meissner solution in the region — L < y < —L + R up to first order in 
r(T) <gC 1 is given by 
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4>(y) — — 4arctan ■ 



H exp 



(L+y) 
Aj 



H s + y/H*-H*' 

H(y) = H 1 (y) = h(y) + H 1 , 



2XjH 



exp 



(L+y) 



H s + y/H* - H 2 + H 2 exp 



2(L+y) 
A, 



j(y) = ji,o(v) = -4H H s + y/Hl - H 2 



[H. + ^H 2 - H 2 


2 

— H 2 exp 


2{L+y)~ 
X J \ 


exp 


(L+y)] 



H s + ^H 2 - H 2 



H 2 exp 



2(L+y) 
Aj 



J{y) = J (y) = My) = — [H-H 1 - h{y)] , 



f(y) EE f (y) = h{y) = 1 



r(T) 



Xj 2 h(y) + -[h(y) + H 1 -HY 



(72) 
(73) 

(74) 



(75) 

(76) 
(77) 



The Meissner solution in the region L — R < y < L can be obtained from (|72j)-(|77j) by means of the substitution 
y — > —y, 4>{y) — * —4>{—y). In the region —R < y < R, the solution is 



(f>(y) = 0, h(y) = 0, j(y) = 0, 
H(y) = H 1 . 



f(y) = i- r -^-[H-H 1 f. 



(78) 
(79) 

(80) 
(81) 



2. The vortex solution 

In the region —R < y < R, the vortex (antivortex) solution satisfies the initial conditions a = </>(0) 
<>n 



= ±2H S [Eq. <Mj] and has the form 



<f>(y) = ±4arctanexp 
h(y) = ±Aj cosh _1 
j(y) = T2cosh^ 2 A 



y_ 



JJ_ 



sinh 



y_ 

Aj 



= 7T, /3 EE 

(82) 
(83) 
(84) 
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The quantities H(y), J(y) and f(y) are given by fl73|), ( ff6| ) and (77|), respectively, with h{y) taken from Note 
that the field induced by a vortex at y — 0, in agreement with (j5C), is 



2 rr 

/i(0) = Xj = 



2 + e 2 ' 



and not h(0) = H s , as in the case of a single junction with thick electrodes. j|] By inserting ( |67j ) and (]82|) into (|54|), 
we obtain the vortex self-energy: 

£„ = 8r(T)-=4=. (85) 
V2 + e 



$ = 7T- 



The vortex flux, according to (b2h, is 



and the lower critical field, by (p5f), is 



2 + e 2 



„ A _ i„, . I^±Z. (86) 

Thus, for e <C 1, the vortex flux <f> <C $0 = ancl the lower critical field (^) is much larger than the corresponding 
field of a single junction with thick electrodes H cW — \ \p^, in agreement with Ref. 0. 

B. A double-junction stack (N — 3) 



In the double-junction case, the nonzero matrix elements (AE) of G(n,m) are 



G(1,1) = G(2,2)= ^±4? , G(l,2) =G(2,1) = ^ . (87) 

(2 + e 2 ) 2 -l (2 + e 2 ) 2 -l 

A 2 A 2 

The corresponding 2x2 matrix G(n, m) (see Appendix A) has two positive eigenvalues: and ^f-, with the lengths 

Aji = —7===, A 2 = 



VT+e 2 V3+^ 
According to (E|), the superheating (penetration ) field is 

H s = Xjl = ^±Z, (89) 

which is smaller than the corresponding single-junction value (|7l|), in agreement with Ref. |l. The application of J2H| ) 
yields the value of the local field in the homogeneous Meissner state: 

^1=^2 = ^. (90) 



1. The Meissner solution 

The Meissner solution in the fields < H < H s = ^ 1 ~^ e ' 2 obeys the symmetry relations (^lj). The substitution of 

My) = My) = <t>{y) (9i) 

into @, using @, yields 

^M) 1 If \ / 00 N 

2- = T2- sll W)' (92) 



dy 2 A 



.71 
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Thus, all the results of the single-j unction case, Eqs. (|7^)-(|8l|), apply if we substitute Aj — > Xj\ (note that Xj < Aji), 
make the identification 

<p( y ) = My) = My), H{y) = H l {y)= H 2 {y), h{ y ) = h l {y)=h 2 {y), 



j(y) = h,o{y) = J2,i(y), J(y) = My) = My), f(y) = My) = My), 

and take the values of H s and Hi from (^) and (|9(]), respectively. Moreover, 

My) = o, 

and 

My) = i- r -^h( y ). 
A J1 



(93) 



(94) 



2. The vortex-plane solution 



In the region —R < y < R, the vortex-plane (antivortex-plane) solution describes two vortices (antivortices) [one 

il(0 

dy 



vortex (antivortex) per I-layer] and satisfies the initial conditions ol\ = 0i(O) = 7r, a 2 = 4> 2 (0) = it, /3i = d ^ >1 ^ — 



±2H S , I3 2 



_ d<fe(0) _ 
— dy 



±2H S [Eq. (E|)]. By (HJ), it obeys the symmetry (91) and Eq. (62), with 



0(y) = ±4arctanexp 



v/i 



Thus, explicit expressions for h\(y) — h 2 (y) = h(y) and ji.o(y) = j 2 ,i(y) = j(y) can be obtained from single- 
junction Eqs. (|83]), d&j), taking account of the substitution Aj — > Aji. The quantities H\(y) — H 2 (y) = H(y), 
My) = My) = J (y) and My) = My) = f(v) are g iven b y ©: @ and @> respectively, with i?i taken from 
( |90| ) . For Ji (y) and /1 (j/) , we have (|9|) and ( pi| ) , respectively. 
The vortex-plane self-energy is 



= 16r(T)Aji = 16r(T) 



VTTe 2 " 1 



(95) 



and the flux is 



$ = 2tt 



1 + e 2 ' 



(96) 



which immediately leads to the lower critical field: 



Hd — —H.i — 



2 vT+e 2 



As can be seen by comparing (|95|) with the single-junction expression flSq) , the energy per vortex in the double-junction 
stack is higher. Finally, the field induced by the vortex plane at y = 0, according to (50), is 



h(0) = hi(0) = h 2 (0) 



e 2 H s 
1 + e 2 
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3. The vortex-antivortex solution 



As can be easily seen, equations (|24|) with N = 3 admit in the region —R<y<R another exact topological 
solution, namely an incoherent vortex-antivortex solution 

Mv) = -My) = Hy), (97) 

where <f>(y) is given by 



4>(y) = 4arctanexp 



JJ_ 

A 2 



With oti — a.2 — 7r and (3\ = —(i 2 = i2A 2 , the vortex-antivortex solution explicitly satisfies the existence criterion 
(Q) . By (JU) , the self-energy of the vortex-antivortex solution is 

E va = 16r(T)A 2 = 16r(T) yJ= , (98) 

V3 + e^ 



which is lower than the self-energy of the vortex-plane solution (|95[) . 

However, the vortex-antivortex solution docs not satisfy integrodifferential equations (j2Fj) with N = 3: These 
equations do not possess the symmetry ( B7| ) for any H ^ 0. The vortex-penetration conditions also cannot be met, 
because /3i/?2 < 0- Moreover, the Gibbs free energy ( p^ ) of the vortex-antivortex pair is always positive with respect 
to the Gibbs free energy of the Meissner state: The flux $ carried by this pair is exactly equal to zero. Thus, the static 
vortex-antivortex solution is thermodynamically unstable, in agreement with the general consideration in section III. 
As shown in Ref. Q, the vortex-antivortex solution can be realized in the dynamic regime, in the presence of an 
external current applied to the central S-layer. 



4- Single-vortex solutions 

Consider a configuration with a single vortex in one I-layer (say, with n = 2) and no vortices in the other. The 
solution representing this configuration does not possess any symmetry. As shown in Appendix C, the self-energy of 
a single vortex, E sv , satisfies the exact inequality E va < E sv < E v , where E v and E va are given by (|95| ) and d9§|), 
respectively. On the other hand, the total flux carried by a single vortex, according to ( |66| ) with N — 3 and I = 2, is 

$ = 7T- 



1 + e 2 ' 

which is exactly half the total flux of the vortex plane (Q). Hence, the Gibbs free energy of a single vortex is positive 
with respect to that of the Meissner state for H < H c \ = — v/1 + £2 . This example clearly illustrates thermodynamic 
instability of single- vortex solutions discussed in section III. 

The two second-order differential equations describing the single-vortex configuration can be reduced p"2| to one 
fourth-order equation 

d 4 d> 2 2 + e 2 . , fd6 2 \ 2 2 + e 2 , d 2 6 2 



■ sin q>2 



dy 4 e 2 V dy I e 2 dy 



d 2 4> 2 
dy 2 



1 - 



(2 + e 2 )sin^ 2 



,2 d 2 <p2 
dy 2 



c 2d s <t>2 



+ t/1 



(2 + e 2 ) sin 02 - e : 



dy 2 



2 + e 2 -1 



■ sin (p 2 



dy 2 



(99) 



The phase difference 4>\{y) can be found without any additional integration from the relations 
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!>i = arcsm 



e)sm0 2 -e — 



7T 7T 

2 - V1 ~ 2 ' 



(2 + e 2 ) sin 02 - e : 



d 2 4> 2 
dy 2 



< 01 < 7T, 



arcsm 



(2 + e 2 ) sin 02 - e 



2 d 2 2 



dy 2 



— 7T, — 7T < 01 < — — . 



The initial conditions for the vortex solution are given by the relations a\ = 0i(O) = 0, a 2 = 02(0) = 7r, /3i 

(2 + e 2 ) 2 -l" 



(100) 



0, /?2 = > and must satisfy the existence criterion ( |&4| ) 



(2 + e 2 ) (/3 2 +/3 2 )-2|/3 1 |/3 2 = 



In view of the condition /?2 > 0, the appropriate solution of (101) for @2 is 

dMO) 



dy 



1 

2 + e 2 



\fh\ + \ 



(2 + e 2 ) 2 - 1 [4 (2 + e 2 ) - e 2 /? 2 ] 



dy 



< 



(101) 



0< < 



2X/2+1 2 



By the symmetry 



d ^i - 1 = 0. The initial condition on 4r^r can be obtained from the relation 

ay z dy A 



d 3 Mv) 
dy 3 



(2 + e ) cos 02 (y) — — cos 0i (y) 



dy 



dy 



that follows directly from (E 



d 3 2 (O) 
dy 3 



(2 + e 2 Y 



1 



[4(2 + e 2 )-e 2 /3 2 ]. 



Thus, we have obtained a complete formulation of the single-vortex problem in terms of the standard initial value 
problem. In agreement with general consideration of section III, the problem admits a family of single- vortex solutions 
parameterized by 

Unfortunately, there are no general methods of analytical integration of nonlinear differential equations of order 
higher than two. However, numerical integration of (^9|) with the above-derived initial conditions should pose no 
problem. Moreover, it is not difficult to obtain asymptotics of the single- vortex solution in the region \y\ ^> Xj%. For 
y <C — Aji, equations (199) , (10C) can be linearized: 



d 4 2 2 (2 + e 2 )d 2 2 | (2 + e 2 ) 2 -! 
dy 4 e 2 



dy 2 



(2 + e 2 ) 2 - e 



0. 



;rf 2 02 

dy 2 



(102) 



(103) 



The solution of (102), (|103J) obeying the vortex asymptotic conditions is straightforward: 

y 



<t>i(y) = Ci(/3i)exp 



A 



.71 



C 2 (/?i)exp 



y_ 

a 2 



fo(y) = Ci(/3i)exp 



Aji 



C 2 (Pi) exp 



y_ 

A 2 



y < -Aji, 



(104) 
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where Ci(/3i),C2(/3i) > are constants with respect to y, parameterized by the initial value \f3\\. The asymptotics 
for y ^> Aji can be obtained from (104) by the use of the symmetry relations (p7h: 



0i (j/) = -Ci(/3i)exp 



V 

A J! 



C 2 (/3i)exp 



y_ 

A 2 



02 (y) =2Tr-Ci(/3i)exp 



y 



C 2 (/3i)exp 



"A 2 



y > Aji, 



(105) 



Expressions (104), (105) illustrate a very important property of single-vortex solutions in multilayer structures: Their 
spatial dependence is characterized, in general, by N — 1 different length scales, which agrees with the conclusions of 
Ref. §. 



C. A 3-junction stack (JV = 4) 



The consideration of all solutions to (|24|) in the case N > 4, including thermodynamically unstable incoherent 
vortex configurations, requires the use of (N — 1) > 3 second-order nonlinear differential equations, or, equivalently, 
of one nonlinear equation of order 2 (N — 1) > 6. Therefore, for the sake of mathematical simplicity, from now on we 
concentrate only on Meissner solutions and topological vortex-plane solutions, both obeying the symmetry (^). (The 
analysis of incoherent vortex configurations for TV > 4 can be done using the algorithm worked out in the previous 
subsection.) 

For N — 4, the Meissner solutions and the vortex-plane solutions satisfy the relation 03 (y) = 0i(y). The relevant 
two equations for 4>i(y) and 2 (y) can be given the form 



rf 4 02 | 2 + £ % m ^ (#2 

dy A e 2 \ dy 



2 + e 2 , d 2 2 



2+e 



sin ( 



d 2 4>2 
dy 2 



(2 + e 2 )cos0 2 ^- e 2 ^ 



(2- 



c 2 d 2 4>2 
dy 2 



1- 



(2 + e 2 )sin0 2 -e 2 



dy 2 



(2 + e 2 ) -2 . 1 2 + e 2 d 2 2 
? Sm<t>2 ~—^ 



(106) 



9i = arcsm — 
2 



(2 + e 2 ) sin 2 



2 d 2 2 



dy 2 



TT TT 

2 - V1 ~ 2 ' 



arcsm — 



(2 + e 2 ) si 



_> . 2 <i 2 02 



sin 02 



dy 2 



TT, ~ < 01 < TT, 



arcsm — 



e I sm 02 



dy 2 



-TT < 01 < --. 



(107) 



Note that Eqs. ( |1 06| ) , (|107| ) have the same mathematical structure as equations for single vortices in the case of a 
double-junction stack (99), (10C). According to (Eq), the superheating (penetration) field for N = 4 is 



H„ = 



(2 + e 2 ) 2 -2 



eVlO + 3e 2 

The local field in the homogeneous Meissner state, according to (E3h, is given by the relations 



Hi = H 3 



(2 + e 2 ) H 
(2 + e 2 ) 2 -2' 



2H 



(2 + e 2 ) z -2 
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1. The Meissner solution for < H <C H s 



In fields < H <§C H s , the general criterion of the existence of Meissner solutions ( fl4| ) takes the form 



a?' 



where a\ = (f>\{—L) -C 1, a| = <j)\{—L) <C 1. Thus, equations (UOq), (107) can be linearized: 



d*^a 2 (2 + e 2 ) d 2 2 (2 + e 2 )' 



1 



dy 2 



-02 = 0, 



(108) 

(109) 
(110) 



The Meissner solution of ( |109| ), (|110|) in the region — L < y < —L + R, obeying the boundary conditions 



My) = -^S 



V2 



(V2 + 1) 



\/2 - y/2 + e : 



■. exp 



(y + L) 



A 



71 



(V2-1) 



\/2 + V2 + e ; 



: exp 



(y + L) 



A 



.72 



b 2 (y) = -eH 



(V2 + 1) 



\/2- V2 + e 2 



exp 



A 



.71 



\/2 + V2 + e 2 



exp 



A 



.72 



where 



Aji 



y/2- \/2 + e 2 



>J2 



V2+72+^ 2 



(111) 



(112) 



We observe that spatial dependence of the solution (111) is characterized by two different Josephson lengths, Aji and 
Aj2, in agreement with the general consideration of section III. Moreover, 4>i(y), 4>2(y) < 0, and \<f> i (y) \ < \4>2(y)\- The 
values a.\ = <f>i(—L) and a 2 = <j>2{—L) meet the existence criterion ( |108| ), as they should. Using (111), we can derive 
explicit expressions for all physical quantities of interest from general formulas of section II. The Meissner solution 
in the region L — R < y < L can be obtained from ( [Hip by means of the substitution <f)i,2{y) — * — 0i, 2 (~2/)- I n the 
region —R < y < i?, we have 4>i 2 (y) = 0, and H n (y) = H n , as usual. 



2. The vortex-plane solution 



a 2 



02(0) = 7T and Pi 



_ d<p 2 (o) 

dy 



The vortex-plane solution to (106), (|107|) in the region —R < y < R obeys the initial conditions u\ 



2H S , fo 



(702(0) 
dy 



l(0) 



2H S [Eq. (fl9|)]. By the symmetry (f46|), we also have 



d < ^°- ) — 0. The initial condition on ^4r is derived from the relation 



dy 2 



dy 3 



d 3 My) 

dy 3 



dy dy 



that follows from (p 



dy 3 



= -2H, 



In this way, we arrive at a complete formulation of the initial value problem for the vortex-plane configuration. In 
contrast to the single-vortex problem considered in the previous subsection, the above-derived initial conditions do 
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not contain any arbitrariness, in full agreement with our general conclusion in section III about the uniqueness of the 
vortex-plane solutions. 

Although the vortex-plane problem admits only numerical integration, the asymptotics in the regions \y\ 3> Xji 
and \y\ <C Aj 2 can be readily obtained. Thus, on the basis of linearized Eqs. (109), ( 1 1C ) , we have: 



1 



C\ exp 



y 



Aji 



C 2 exp 



V2 



h(y) = Ci exp 



Aji 



+ C 2 exp 



V/2 



y < -Aji; 



0i(l/) = 2tt 



V2 



Ci exp 



y 



C 2 exp 



Aj2 



6 2 (y) = 27r - Ci exp 



Aj! 



C 2 exp 



y 

Aj 2 



2/ > A 



.71. 



where C\ and C 2 are positive constants. The solution in the region \y\ <C Aj 2 is represented by Taylor series expansions: 



My) = n + 2H s y 



:k 2 



-y 3 



Mv) = n + 2H s y- -H s y 3 



Finally, the total flux carried by the vortex plane, by (52), is 

_e 2 (l0 + 3e 2 ) 



$ = 7T 



(2 + e 2 ) 2 -2' 



and the exact value of the induced field at y = 0, by (pQ), is 



(See Fig. 2.) 



e 2 3 + e 2 )H S e 2 4 + e 2 H s 

h^O) = h 3 (0) = — K - J—l h 2 (0) = — i 

W W (2 + e 2 ) 2 -2 (2 + e 2 ) 2 -2 



D. A 4-junction stack (iV = 5) 



For N ~ 5, the Meissner solutions and the vortex-plane solutions satisfy the relation cj>4 (y) 
6 2 (y). The equations for 4>i(y) and </> 2 (y) reduce to the form 



My) and ^3(2/) 



1 + 6 2 



dy 4 



■ sin <p 2 



1 + e 2 



■ COS 2 ■ 



dy 2 



sin 2 



d> 2 



(l + e2)cos ^ 2 ^_ e 



dy 3 



(1 + e 2 ) sin0 2 



' ~~d!F 



(1 + e 2 ) sin ^2 



, d 2 2 
dy 2 



(2 + e 2 ) (l + e 2 )-l 



e 2 d 2 2 



dy 2 



(113) 
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!>i = arcsm 



e 2 ) sin d>2 - e 2 



d 2 02 



dy 2 



7T 7T 

-- < d>! < -, 
2 - V1 - 2 ' 



arcsm 



'l + e 2 ) sin 02 - e 2 



rf 2 02 
dy 2 



7T 

2 < 



< 7T, 



arcsm 



. 2 d 2 ^ 2 



l + e z sin 02 -e 



dy 2 



7T, — 7T < 01 < — — . 



(114) 



Mathematical structure of Eqs. (113), (114) is analogous to that of Eqs. (106), ( |107| ) of the 3-junction stack, which 
allows us to skip some details in what follows. According to (45), the superheating (penetration) field for N — 5 is 



_ V2[(2 + 6 2 )(l + e 2 )-l] 
eV5 + 2e 2 

The local field in the homogeneous Meissner state, according to (p3|), is 

(l + e 2 )i/ 



-Hi = i?4 



(2 + e 2 )(l + e 2 )-l : 



Ho — H'i 



H 



(2 + e 2 )(l + e 2 )-l 



1. The Meissner solution for < H <C H s 
For the fields < H <C H s , the Meissner solution in the region —L<y<—L + R has the form 



My) = - 



eH 
2\7I 



5 - 1) (3 + V5) 



3-V5 



exp 



(y + L) 



A 



.71 



(V5 + 1) (3 -a/5) 



3+V5 . ,2 
2 ' 



exp 



A 



.72 



where 



(3 + V5) 

3-V5 , ,2 
2 ' c 



Aji = 



: exp 



Aji 



3-\/5 I ,2 
2 T c 



=, A 



.72 



(3 -a/5) 

3+yg , 2 
2 "t" c 



: exp 



A 



72 



3+75 , ,2 
2 ' c 



(115) 



By comparison with Aji, Aj2 of the 3-junction case, Eq. ( |112| ), the Josephson lengths given by (115) are larger, 
which illustrates a general tendency: Josephson lengths increase with increasing N. As in the 3-junction case, 
0i (y) , 4>2 (y) < 0, and |0i(t/)| < 1 02 (?/) I - The values a% = 0i(— L) and a 2 = 02 satisfy the existence criterion 
(pi), as expected. 



2. T/ie vortex-plane solution 



The vortex-plane solution to (113), (114) in the region —R<y<Ris characterized by the initial conditions 
(0) = 7T, a 2 = 02 (0) = 7T and ft = = 2ff s , /3 2 = = 2 i7 s [Eq. ©]. In addition, we have 

, r . — ■ ,„„, , ; - -2H S , as in the 3-junction case. 
The asymptotics of the vortex-plane solution in the region \y\ 3> Xji are 



d 2 4>2(0) _ q and d 3 <fe(0) 



i (2/) 



\/5-l)Ci 



exp 



\/5 + 1] C 2 exp 



V2 
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h(y) = Ci exp 



Aji 



+ C 2 exp 



V72 



y < -Aji; 



Si(»)=2tt-- 



5 — 1 1 C\ exp 



Aji 



1 ) C 2 exp 



y 

Aj2 



02 (y) = 27r - Ci exp 
where Ci, C 2 > 0. In the region \y\ Aj 2 , we have: 



y 

A./i 



C 2 exp 



y 

Aj 2 



=n + 2H s y 



[l + e 2 )H s 3 
" V 



y > A 



02 (y) = TT + 2H s y- -H s y 3 



The total flux carried by the vortex plane is 



$ = 2tt 



(5 + 2e 2 ) 



(2 + e 2 )(l + e 2 ) -1' 



and the exact value of the induced field at y = is 

e 2 (2 + e 2 ) if, 



&i(o) = MO) 



(2 + e 2 )(l + e 2 )-l' 



h 2 (0) = h 3 (0) 



e 2 (3 + e 2 ) H s 
(2 + e 2 )(l + e 2 ) -1" 



(See Fig. 2.) 



E. The layered-superconductor limit (TV - 1 > 2 

The limit TV— 1 > 2 [i] ( \-] stands for the integer part of -) corresponds to the situation when a periodic thin-layer 
S/I structure can be regarded as a "layered superconductor" rather than merely a (TV — l)-Josephson-junction stack. 
Indeed, in this limit the superheating (penetration) field, Eq. (Eq), becomes 



I 



e(TV-l) 



and for (TV — I ) — > 00 tends to the limiting value of an infinite layered superconductor H soo — 1. (We remind 
that the lower critical field of an infinite layered superconductor, according to Refs. [QJ|] , is H cloo = -■) The total 
flux carried by a vortex plane, by (p2|), in the considered limit is 



$ = tt(TV - I) 



2a/1 



e(TV-l) 



which for (TV — 1) — ► 00 tends to the limiting value of one flux quantum $0 — t per vortex, as expected for an 
infinite layered superconductor. Moreover, according to (^3|), for I-layers whose index n satisfies the condition 
[7] < n < TV - 1 - [i] , we have 

\H n \ = \H\ (^ + M JV - n )«|ff|. 

In other words, for \H\ < H soo — 1, the inside junctions exhibit the complete Meissner effect in the region —R < y < R 
. Thus, the integer [i] determines the number of junctions near the boundaries x = and x = TV — 1 that "feel" the 
influence of the superconductor/ vacuum interfaces. Far from the boundaries, i.e. in the region [i] <C x <C TV — 1 — [i] , 
one can apply all the results of the theory of infinite layered superconductors. p]|| 
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V. DISCUSSION 



Within the framework of standard methods of the theory of ordinary differential equations, we have obtained 
a complete mathematical description of Josephson vortices and of the Meissner effect in periodic thin-layer S/I 
structures. A rigorous examination of the properties of Eqs. (§3) allowed us to establish the general existence 
criterion (33), which formed a solid basis for our subsequent physical and mathematical conclusions. One of the most 
striking physical consequences is that the derivation of the exact expression for the vortex-penetration field H s , Eq. 
([45|), did not require any explicit solution of (p4j). Although explicit analytical solutions to Eqs. (|^) proved to be 
possible only in a limited number of cases discussed in section IV, numerical integration of these equations should 
pose no problem owing to the algorithm worked out in the paper. 

All the three types of localized solutions obtained in the paper possess a number of interesting physical and 
mathematical properties. For example, the Meissner solutions are characterized by several different Josephson lengths 
^Ji (y lengths for even TV, and ^-=1 lengths for odd N). Unfortunately, this important fact was not noticed in previous 
publications. Jl4j We think that our result may prove to be useful in view of the current experimental efforts Jl6|,|l7j to 
verify the interlayer tunneling model of high-T c superconductivity ]l8|-|2"2|] by measuring the c-axis penetration depth. 
[The penetration of the parallel magnetic field with a distribution of length scales has been recently observed [^3| in 
the organic layered superconductor k-(BEDT-TTF) 2 Cu(NCS) 2 .] 

Mathematically, Josephson vortices, represented by both the vortex-plane solutions and incoherent vortex solutions, 
are static sine-Gordon- type solitons. They satisfy the standard |ll| asymptotic boundary conditions, Eqs. (ff7|), (|4g|), 
and Eqs. (|58|)-(|6C). The expression for their self-energy (|54|) is a direct generalization of the well-known expression for a 
single junction. |8|| The thermodynamically stable vortex-plane solutions demonstrate a profound difference between 
Josephson-vortex formation in weakly-coupled multilayer structures and Abrikosov-vortex formation in continuum 
type-II superconductors (isotropic or not). The proof of the existence of such solutions in the general case of (N — 1)- 
junction stacks establishes relationship to the well-know "coherent mode" (alias the "in-phase" mode) in a double- 
junction stack [|] |7| and the recently obtained Jl],|| vortex-plane solutions in infinite (N = 00) layered superconductors. 
However, in contrast to the latter two cases, the vortex-plane solutions for 4 < N < 00 are characterized by several 
Xji, as the Meissner solutions. Moreover, in contrast to the case N = 00, the intralayer currents, J n , in the presence 
of a vortex plane are not equal to zero, with the exception of J n-i in the case of odd N . 

The single- vortex solutions are not uniquely determined by the asymptotic boundary conditions (|58|)-(|60|), as follows 
from the existence criterion (|64]). Their spatial dependence is characterized, in general, by N — 1 length scales. In 
contrast to the vortex planes, isolated Josephson vortices cannot penetrate the periodic S/I structure at any \H\ 7^ 
and do not satisfy the original integrodifferential equations (||) , ( |20| ) minimizing the Gibbs free-energy functional (Q) . 
Although the self-energy of incoherent vortex solutions may be lower than the self-energy of a vortex plane, their 
Gibbs free energy is positive with respect to the Gibbs free energy of the Meissner state for \H\ < H c \ {H c i is the lower 
critical field for the vortex-plane solution), as is illustrated by our analysis of a double-junction stack in section IV. 
However, incoherent vortex solutions can be realized in the dynamic regime. In layered superconductors, isolated 
vortices may also emerge as metastable topological entities owing to pinning by extended defects. We emphasize that 
our general conclusions about the form of single-vortex configurations stand in full agreement with numerical results 
ofRef. {§. 

Finally, our exact results clearly show that the phase-difference equations ( p4[ ) do not admit any solutions in the 
form of Abrikosov-type vortices, suggested in Refs. As we have already pointed out, the hypothesis 

of Abrikosov-type solutions in infinite layered superconductors leads to incorrect estimates of the value of the lower 
critical field i? c ioo- Unfortunately, it seems that some other theoretical predictions, based on the assumption of the 
existence of Abrikosov-type solutions, should also be revised. 



APPENDIX A: THE SOLUTION OF THE FINITE DIFFERENCE EQUATION FOR H N {Y) 

Equations (^o|) can be regarded as a nonhomogeneous finite difference equation for H n (y) with respect to the layer 
index n, subject to boundary conditions (||). According to general theory of such equations, 0] its solution for e < f 
can be represented in the form 

H n (y) = h n {y) + H n , (Al) 

where 
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JV-1 



is the particular solution of 



and 



K{y) = y X G ( n > '' 



711 — 1 



dy 



satisfying the boundary conditions 
h o{y) = /iiv(y) = 0, 



H (n~ n + n~ N+n - n n - fi N - n ) 



(A2) 



(A3) 



(A4) 



M = l + |-^1+^<1 



(A5) 



is the solution of the homogeneous form of ( J2C)| ) (with the zero right-hand side) meeting the boundary conditions 

H = H N = H. (A6) 

The quantities G(n, m) in ( [A2] ) are matrix elements of (A*" + 1) x [N + 1) matrix Green's function G(0 < n, rn < N). 
They obey the nonhomogeneous finite difference equation 



G{n + 1, m) - (2 + e 2 ) G{n, m) + G(n - 1, m) = -S n 
(fin,m is the Kronecker index) with the boundary conditions 

G(0,m) = G{N,m) = 0. 

The explicit form of G(n, m) is 



G(n, m) 



2eJl 



4 



..n { ..m — N ,.N—m\ i ..N—n ( ,.— m ..Wt 

|n-m| _ /; \P M J t l/t M J 



The following properties of G(n, m) can be easily verified using ( |A7| ) and (|A9|): 

G(n,m) = G(m, n), 



(A7) 



(A8) 



(A9) 



(A10) 



G(n, N — m) — G(N - n, m), 



G(n, m) > for any 1 < n,m < N — X, 



(All) 
(A12) 



JV-l 



J2 G(n, m) = 4 [1 - G(n, 1) - G(n, JV - 1)] 



1< n < JV- 1, 



(A13) 



JV-l JV-1 



AT- 1 



AT 



(A14) 



Note that matrix Green's function for an infinite layered superconductor, (n,m), is determined by the matrix 
elements 



25 



p\n-m\ 

Goo(n,m) — . — oo <n,m < +00, (A15) 



that satisfy the summation rule 

N-l 1 

^2 G co (?i,m) = — . 

m— 1 

Consider now a (N — 1) x (JV — 1) matrix G(n, m), whose matrix elements are given by the right-hand side of 
(A9) with 1 < n, m < N — 1. Of special physical importance are positive eigenvalues of G(n, m): They determine 
characteristic length scales of localized solutions to (|4j). Indeed, in the asymptotic region where = o(l) and 

4>n{y) = o(l), equations ( p7| ) can be linearized with the result 



^ 1 G (n,m)^# = i0 n (y), l<n<iV-l. 



The substitution 4> n (y) oc exp [±f ] yields 



at-i a2 

2 G(n, m)4> m {y) = -^Mv)> 1 < n < N - 1, 

m— 1 

which is exactly an eigenvalue equation for G(n, m), with p- being a positive eigenvalue. 

APPENDIX B: PROOF OF THE LEMMA 

By introducing new functions 

ipi(v) = (t>i{y),^2{y) = <f>2(y), ■ ■ .,ip N -i(y) = <f>N-i{y), 

i 1 s d4>i{y) d<j) 2 (y) . . d4> N - X {y) 

tpN{y) = —, — ,tpN+i{y) = —\ — , • • • , ip2N-2{y) = — -, , (Bi) 

ay ay dy 

we convert (^4|) into an equivalent normal system of 2N — 2 first-order equations 

^™= J P i (Vi,V2,...,^2iV-2), l<*<2AT-2, (B2) 
ay 

Fi (ipi,ip2, ■ ■ -,^2^-2) = ^i+JV-i, 1 < i < N - 1, 

F;v (lpl,1p2, ■ ■ ■ ,1p2N-2) = t( 2 + f2 ) Sm ^l - Sin ^2] , 

Fi (iI>i,i/j2,---,i/J2N-2) = ^ [(2 + e 2 ) sin^i -sin^-i -sin^+i] , N + l<i<N — 3, 

F-iN-2 (lpl,1p2, ■ ■ -,4>2N-2) = "2 [( 2 + £ 2 ) Sm-0AT_i - Sm^iV-2] , 

subject to initial conditions 

Mvo) = oti, 1 < * < iV — 1, 



2G 



ipi(y ) = fo-N+i, N<i<2N-2. 



(B3) 



To prove the statement of the Lemma, it is sufficient to observe that all Fi (ipi , ip2 , ■ ■ ■ 
functions of their arguments for y E (— oo,+oo) and ipk € (— oo,+oo) (1 < k < 2N — 2). 
derivatives with respect to tpk satisfy the relation 



8Fi (lpl,1p2, ■ ■ - ,1p2N-2) 



dipt 



< 



,ip2N-i) are continuous 
Moreover, their partial 



(B4) 



for y E (— oo,+oo) and ipk E (—00, +00) (1 < i, k < 2N — 2). Thus, the Lipschitz conditions with respect to 
ipk are met for y E (—00, +00) and ipk € (— 00, +00) (1 < k < 2A*" — 2), which immediately guarantees the 
existence and uniqueness of a solution to (B2), satisfying arbitrary initial conditions (B3), in an arbitrary interval 
I = [£1,^2] such that yo E I. Continuous dependence of the solution on initial data is a result of continuous 
dependence of Fi (ip%,ip2, ■ ■ ■ ,ip2N-2) on their arguments and of the condition (B4). Infinite differentiability of the 
solution automatically follows from infinite differentiability of Fi (ip\, ip2, ■ ■ ■ , 1P2N-2) with respect to their arguments. 



APPENDIX C: THE UPPER AND THE LOWER BOUNDS FOR THE SELF-ENERGY OF 
SINGLE- VORTEX SOLUTIONS IN A DOUBLE-JUNCTION STACK 

To determine the upper and the lower bounds of the self-energy of single- vortex solutions in a double-junction stack, 
we must find the extrema of the right-hand side of (B4J) with N — 3 under the normalization condition 



2 

E 

n=l 



dy 



d(j>n{y) 

dy 



= const < 00. 



This leads to the variational principle 



S<f>n(y) 



2 2 
k—l m—l 



\ - \ - / (lot,- {ll ) llO„,(ll) 



du 



du 



2 

fe=i 



j du 


~d4> k (u)~ 




du 



where A 2 is a Lagrange multiplier. Performing the variation with the use of the boundary conditions ^"^y^ 
we arrive at the eigenvalue problem for the 2x2 matrix G(n, m), determined by the matrix elements (| 

d 2 (f) m {y) _ A 2 d 2 4> n (y) 



2 

m=l y 



J 2 



y 



The only two eigenvalues of G(n,m) are and ^f, with Aji and A2 given by (pq). The larger eigenvalue, 
corresponds to the vortex-plane solution with <j>i — (f>2 and determines the upper bound ( |95| ) for the self-energy. The 

A 2 

smaller eigenvalue, -p-, corresponds to the vortex-antivortex solution with (pi — —(p2 and determines the lower bound 
( |98| ) for the self-energy. Thus, for the self-energy of the single- vortex solutions, E sv , we necessarily have 



Fya F sv *C E v . 
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FIGURE CAPTIONS 
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1. Fig. 1. The geometry of the problem (schematically). Here N = 12; Xj^ << R < L/2, A Joo = 8Trejop; H > 0. 

2. Fig. 2. The distribution of the self-induced field of a vortex plane h(x,y) at y = 0, for e = 0.5: a) a 3-junction 
stack, H s — 1.849; b) a 4-junction stack, _ff s = 1.662. 
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